program powermardiau
use moments3
use matrixfuns
use datahadler
use linear_operators
use imhoff
implicit none
integer, parameter:: n=3,p=10,num=100,T=5000,pth=9,s=1000
double precision,parameter:: si=.999d0
double precision wmat(p+1,p+1),eta(num),power(num,3),pr,dchiin&
&,b(3),iwmat(p+1,p+1),matm(p+1,1),matv(p+1,p+1)
double precision sigma(n,n),MD(p+pth+1),VD(p+pth+1,p+pth+1)&
&,dmut(pth,n),dvsig(pth,n*n),dmat(n,n),matA(pth,pth),matB(p+1,pth),matD(p+1,pth+p+1)&
&,crit,critk,critsk,meank,vark,dnordf
integer i,j,k,ifault,irank



forall(i=1:pth,j=1:n)
dmut(i,j)=0.0d0
end forall
forall(i=1:pth,j=1:n*n)
dvsig(i,j)=0.0d0
end forall
forall(i=1:n)
	dmut(i,i)=1.0d0
end forall
forall(i=1:n,j=1:n)
	dmat(i,j)=0.0d0
end forall
k=0
do j=1,n,1
	do i=1,n,1
		k=k+1
		dmat(j,i)=1.0d0
		dmat(i,j)=1.0d0
		dvsig(n+1:pth,k)=vech(dmat)
		dmat(j,i)=0.0d0
		dmat(i,j)=0.0d0
	end do
end do


eta=linspace(.00001d0,.02d0,num)
b=vassign(n,0.75d0)
!sigma(1,:)=(/1.0d0,0.1d0,0.2d0/)
!sigma(2,:)=(/0.1d0,1.0d0,0.3d0/)
!sigma(3,:)=(/0.2d0,0.3d0,1.d0/)
sigma=eye(3)


!sigma=eye(3)
crit=dchiin(.95d0,dble(p+1))
critk=dchiin(.95d0,dble(1))
critsk=dchiin(.95d0,dble(p))


call var0u(wmat,dmut,dvsig,pth,sigma,b)
iwmat=.i.wmat

do i=1,num,1

	call mardiamvu(MD,VD,matA,matB,b,eta(i),si,sigma,dmut,dvsig,pth)
	matD(:,pth+1:pth+p+1)=eye(p+1)
	matD(:,1:pth)=matmul(matB,.i.matA)

	matm(:,1)=dsqrt(dble(T))*matmul(matD,MD)
	matv=(matD.x.VD).xt.matD

	!Joint
	call imhofprob(pr,ifault,matm,matv,iwmat,crit)
	power(i,1)=1.0d0-pr

if (ifault .ne. 0) then
	print*, "joint",ifault
end if
!Skewness
	call imhofprob(pr,ifault,matm(2:p+1,:),matv(2:p+1,2:p+1),.i.wmat(2:p+1,2:p+1),critsk)
	power(i,2)=1.0d0-pr
if (ifault .ne. 0) then
	print*, "skewness",ifault
end if
!Kurtosis
meank=matm(1,1)/dsqrt(wmat(1,1))
vark=matv(1,1)/wmat(1,1)

power(i,3)=1.0d0-dnordf((dsqrt(critk)-meank)/dsqrt(vark))+dnordf((-dsqrt(critk)-meank)/dsqrt(vark))

	print*, i,real(power(i,:))

end do
call exportvec(eta,"ema.txt",1)
call exportvec(vec2(power),"powrma3a.txt",1)


end program powermardiau